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We analyze a specific class of random systems that are driven by a symmetric Levy stable noise. 
In view of the Levy noise sensitivity to the confining "potential landscape" where jumps take place 
(in other words, to environmental inhomogeneities) , the pertinent random motion asymptotically 
sets down at the Boltzmann-type equilibrium, represented by a probability density function (pdf) 
p*(x) ~ exp[— $(a;)]. Since there is no Langevin representation of the dynamics in question, our 
main goal here is to establish the appropriate path-wise description of the underlying jump-type 
process and next infer the p(x,t) dynamics directly from the random paths statistics. A priori given 
data are jump transition rates entering the master equation for p(x,t) and its target pdf p*{x). We 
use numerical methods and construct a suitable modification of the Gillespie algorithm, originally 
invented in the chemical kinetics context. The generated sample trajectories show up a qualitative 
typicality, e.g. they display structural features of jumping paths (predominance of small vs large 
jumps) specific to particular stability indices p € (0,2). 

PACS numbers: 05.40.Jc, 02.50.Ey, 05.20.-y, 05.10.Gg 



I. INTRODUCTION 

Many random processes in real physical systems admit 
a simplified description based on stochastic differential 
equations. In such case there is a routine passage proce- 
dure from microscopic random variables to macroscopic 
(statistical ensemble) data. The latter are encoded in the 
time evolution of an associated probability density func- 
tion (pdf) which is a solution of a deterministic transport 
equation. A paradigm example is the so-called Langevin 
modeling of diffusion-type and jump-type processes. The 
presumed microscopic model of the dynamics in exter- 
nal force fields is provided by the Langevin (stochastic) 
equation whose direct consequence is the Fokker-Planck 
equation, [l[ and @. We note that in case of jump-type 
processes the familiar Laplacian (Wiener noise generator) 
needs to be replaced by a suitable pseudo-differential op- 
erator (fractional Laplacian, in case of a symmetric Levy- 
stable noise). 

We pay a particular attention to jump-type processes 
which arc omnipresent in Nature (see |3[ and references 
therein). Their characterization is primarily provided by 
jump transition rates between different states of the sys- 
tem under consideration. However our major focus is on 
a specific class of random systems which are plainly in- 
compatible with a straightforward Langevin modeling of 
jump-type processes and, as such, are seldom addressed 
in the literature. 

To this end we depart from the concept, coined in an 
isolated publication [j] , of Levy flights-driven models of 
disorder that, while at equilibrium, do obey detailed bal- 
ance. The corresponding research line has been effec- 
tively initiated in Refs. It has next been expanded 
in various directions, with a special emphasis on so-called 
Levy-Schrodingcr semigroup reformulation of the orig inal 
probability density function (pdf) dynamics, @, Il5l Il6j 
and i-El, c.f. also [TMH. We note in passing that 
the familiar Fokker-Planck equation can be equally well 



formulated in terms of the Schrodinger semigroup and 
this property is universally valid in the standard theory 
of Brownian motion, [l|, [Tj] . Its generalization to Levy 
flights is neither immediate nor obvious. It is often con- 
sidered in the prohibitive vein following [l?], • 

In fact, in relation to Levy flights, a novel fractional 
generalization of the Fokker-Planck equation has been 
introduced in Refs. to handle systems that are 

randomized by symmetric Levy-stable drivers. In this 
case, contrary to the popular lore about properties of 
(Langevin- based) Levy processes c.f. Refs. [l7(-[lll an< ^ 
|20|, the pertinent random systems are allowed to relax 
to (thermal) equilibrium states of a standard Boltzmann- 
Gibbs form. 

The underlying jump-type processes, in the station- 
ary (equilibrium) regime, respect the principle of detailed 
balance by construction Their distinctive feature, if 
compared with the standard Langevin modeling of Levy 
flights, is that they have a built-in response not to exter- 
nal forces but rather to external force potentials. These 
potentials are interpreted to form confining "potential 
landscapes" that are specific to the environment. Levy 
jump-type processes appear to be particularly sensitive 
to environmental inhomogeneities, [a. Il2|. 

Levy flights are pure jump (jump-type) processes. 
Therefore, it seems useful to indicate that various model 
realizations of standard jump processes (jump size is 
bounded from below and above) can be thermalized by 
means of a specific scenario of an energy exchange with 
the thermostat. It is based on the principle of detailed 
balance. We have discussed this issue in some detail be- 
fore [l3j along with an extension of this conceptual frame- 
work to Levy-stable processes. Not to reproduce easily 
available arguments of past publications, we shall be very 
rudimentary in our motivations. 

We quantify a probability density evolution, compati- 
ble with a jump- type process on R (this limitation may in 
principle be lifted in favor of R n ), in terms of the master 
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equation: 



d t p(x) 



[w <j> (x\y)p(y)-w cj> (y\x)p(x)]dy, (1) 



L<\x — y\<£2 



where £1 and £2 are, respectively, the lower and upper 
bounds of jump size and 



w<t>{x\v) 



C, 



cxp[($(y) - $(z))/2] 



\x - y\ 1+fJ - 
r(l + p) sin(7r/i/2) 



(2) 



is the jump transition rate from y to x. We stress that 
w^(x\y) is a non-symmetric function of x and y. 

An implicit Boltzmann-type weighting involves a 
square root of a target pdf p*(x) ~ exp[— and ac- 
counts for the a priori prescribed "potential landscape" 
whose confining features affect the jump-type pro- 
cess. What matters is a relative impact of a confinement 
strength of <I>(x) (level of attraction, see Ref. Q) upon 
jumps of the size \x — y\, both at the point of origin y 
and that of destination x. In principle, $(x) may be 
an arbitrary function that secures a normaliza- 
tion of exp(— In this case, the resultant pdf is 
a stationary solution of the transport equation ([1]) with 
unbounded jump length, e.g. £\ — > and £2 — > 00. 

We note that the presence of lower and upper bounds 
of the jump size £1,2, that are necessary for an implemen- 
tation of numerical algorithms, enforces a truncation of 
the jump-type process (without any cutoffs) to a stan- 
dard jump process. The transition rates of the latter, 
however, are ruled by Levy measures of symmetric Levy 
stable noises with p € (0, 2). A lower bound for the jump 
size is usually removed while evaluating the correspond- 
ing integrals in the sense of their Cauchy principal values. 
An upper bound is less innocent and its effects need to 
be controlled by long tailed pdfs which stands for a dis- 
tinctive feature of Levy flights, see a discussion of Levy 
stable limits of step processes in Ref. Q. There is also 
pertinent discussion of a long time behavior of (uncon- 
fined, e.g. free) truncated Levy flights in Ref. |24j . 

In contrast to procedures based on the Langevin mod- 
eling of Levy flights in external force fields, [1, [H, [l9j . 
there is no known path-wise approach underlying the 
transport equation ([1]). With no direct access to sam- 
ple trajectories of the stochastic process in question, a 
method must be devised to generate random paths di- 
rectly from jump transition rates @. The additional 
requirement here is that we set a priori a " potential land- 
scape" $(x) for a chosen jump- type (symmetric Levy sta- 
ble) noise driver. 

The outline of the paper is as follows. First we describe 
our modification of the Gillespie algorithm which entails 
a numerical generation of random paths for the dynam- 
ics determined by Eqs ([TJ) and ©. Next the statistics 
of random paths is addressed and various accumulated 



data are analyzed with a focus on inherent compatibility 
issues. 

We analyze generic (Cauchy, quadratic Cauchy) and 
non-generic (Gaussian and locally periodic) examples of 
target pdfs for the jumping dynamics. Random paths are 
generated in conjunction with representative Levy stable 
drivers, like e.g. those indexed by p = 1/2, 1, 3/2. Their 
qualitative typicality is emphasized. 

Statistical data, acquired from our modification of 
Gillespie algorithm, have been employed to generate 
the dynamical patterns of behavior p{x,t) — > p*(x), to 
demonstrate the compatibility of the transport (master) 
equation ([TJ , ^) and its underlying path- wise representa- 
tion. Both coming from the predefined knowledge of the 
target pdf and non-symmetric (biased) jump transition 
rates. 



II. RANDOM PATHS: MODIFIED GILLESPIE 
ALGORITHM. 

Here we adopt [25[ (and properly adjust to handle 
Levy flights) basic tenets of so-called Gillespie's algo- 
rithm |2l[ |22j. Originally, this algorithm had been de- 
vised to simulate random properties of coupled chemical 
reactions. The advantage of the algorithm is that it per- 
mits to generate random trajectories of the correspond- 
ing stochastic process directly from its (jump) transition 
rates, with no need for any stochastic differential equa- 
tion and/or its explicit solution. We emphasize that this 
feature of Gillespie's algorithm is vitally important, since 
Langevin modeling is not operational in our framework. 

We rewrite Eq. (UJ in the form [x — y = z) 



d t p{x) = 



El<\z\<62 



-Wj>(z + x\x)p(x) 



Wj,(x\z + x)p(z + x) — 



dz. 



(3) 



To construct a reliable path generating algorithm consis- 
tent with Eq. (J3J we first note that chemical reaction 
channels in the original Gillespie's algorithm may be re- 
interpreted as jumps from one spatial point to another, 
like transition channels in the spatial jump process. An 
obvious provision is that the set of possible chemical re- 
action channels is finite (and generically low), while we 
are interested in all admissible jumps from a chosen point 
of origin xq to any of [xq — £2, xq — £1] U [xq + £1, xq + £2]- 
It is clear that such jumps form an infinite continuous 
set. With a genuine computer simulation in mind, we 
must respect standard numerical assistance limitations. 
Surely we cannot admit all conceivable jump sizes. As 
well, the number of destination points, even if potentially 
enormous, must remain finite for any fixed point of ori- 
gin. 

Our modified version of the Gillespie's algorithm, ap- 
propriate for handling of spatial jumps is as follows [26[: 
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(i) Set time t = and the point of origin x — xq. 



and W(xo) = W 1 (x ) + W 2 (x ). 



(ii) Create the set of all admissible jumps from xq to 
xq + z that is compatible with the transition rate 
w^{z + x \x ). 



(iii) Evaluate 



Wi(x ) 
W 2 (x ) 



W4,(z + x \xo)dz : 



w^(z + x \x )dz 



(4) 



(iv) Using a random number generator draw p £ [0, 1] 
from a uniform distribution. 



(v) Using above p and identities 



/ w^(z + x \x )dz = pW(x ), p <Wi{x )/W(x ); 

— S2 

6 

Wi(x ) + J w^iz + x \x )dz =pW(x ), p>Wi(x )/W(x ), 

El 

I 



(5) 



find b corresponding to the "transition channel" 

Xq -> b. 

(vi) Draw a new number q € (0, 1) from a uniform dis- 
tribution. 

(vii) Reset time label t = t + At where At = 
-\nq/W(x ). 

(viii) Reset xq to a new value Xq + b. 

(ix) Return to step (ii) and repeat the procedure anew. 

Comment 1: The original Gillespie algorithm em- 
ploys a discrete label v (with a finite range) enumerating 
possible chemical reactions channels. To identify a chan- 
nel, one must look for estimates of a double inequality 
(see Eq. (21b) of Ref. H3) 



E 



M 



a Q 



(6) 



i/=i 



where r 2 is a random number, M indicates a total num- 
ber of chemical reaction channels and a v dt stands for a 
probability that the f-th reaction would actually take 
place in the interval (t,t + dt). To adjust this recipe 
to our settings, we need to enumerate the infinite num- 
ber of (infinitely close) channels. This corresponds to 
passing from summation to integration in Eq. (|6]). As 
it has been pointed out above, such situation corre- 
sponds to possible transitions from xq into an interval 
[xo — £ 2 ,xq — £i] U [xq + si,xq + 62]- As the Lebesgue 
measure of a point equals zero, we can replace inequal- 
ities by identities in ©, see step (v) of the above algo- 
rithm. Formally, we can say that although the number 
of jumps destinations is finite, their number is so large 



that it is consistent to approximate it by a dense subset 
of intervals [xq — E2,xq — e{\ U [xq + Ei,xq + £2]- This 
justifies a replacement of finite sum by an integral over 
corresponding interval. The following jump size bounds 
(integration boundaries) were adopted in the numerical 
procedure: S\ = 0.001 and £2 = 1- 



III. STATISTICS OF RANDOM PATHS: PDF 
TIME EVOLUTION AND COMPATIBILITY 
ISSUES. 

Our main task in the present section is to select, 
to some extend generic, transition rates for jump-type 
processes that will prove to be amenable to the out- 
lined random path generation procedure. Once suit- 
able path ensemble data are collected, we shall verify 
whether statistical (ensemble) features of generated ran- 
dom trajectories are compatible with the master equation 
(JXJ) . That includes a control of an asymptotic behavior 
p(x,t) — > p*(x) when t — > 00. 



A. Harmonic confinement (Gaussian target) 

Let us consider an asymptotic invariant (target) pdf in 
the Gaussian form: 



P*0) = -j= e x ■ 
The corresponding /i-family of transition rates reads 

e -z 2 /2+xz 



w^(z + x\x) = C M - 



i+it 



(7) 



(8) 



4 



According to the step (iii) of the simulation algorithm, 
we must evaluate integrals with transition rates ([5J in 
the intervals [— £2,— £1] and [ei,£2]- To execute the step 
(iv) of the algorithm, we employ the Mersenne- Twister 
random number generator, [23( . We find b by a numerical 
solution of the transcendental equation ([5]) in conformity 
with step (v) of the algorithm. The C-codes for trajectory 
generating algorithm, [26j . were finally employed to get 
the trajectory statistics data for three specific choices of 
Levy drivers, namely /i = 0.5, 1, 1.5. 

Comment 2: For small z, transition rates vary 
rapidly so that adaptive numerical integration algorithms 
become very time-consuming. To speed up the calcula- 
tions, we propose a more efficient procedure (than adap- 
tive numerical integration algorithms with huge number 
of subdivisions for small z), that amounts to separate 
integrations for "large" and "small" z subintervals. At 
small z we can expand the corresponding integrand in 
Taylor series and truncate it at, say, quadratic term. 
Generally speaking, the number of terms to be left de- 
pends on the accuracy which should be retained during 
the integration. At "large" z the transition rates vary 
gradually so that standard adaptive numerical integra- 
tion works fine. 

As an example, we consider the [e\ 1 £2] integration with 
fjt=l. We set £12 = 0.05, so that for z € [£1, £12] we get 



paths, we have not found significant qualitative differ- 
ences to justify a presentation of data for 100 000, 200 
000, 250 000 and more trajectories. The relaxation time 
rate dependence on /i is clearly visible as well. It suffices 
to analyze differences between three curves for t = 0.2 
and/or t = 1. We observe a conspicuous lowering of their 
maxima with the growth of fi (take care of different scales 
on the vertical axes on Fig. 2 panels). The simulated pdfs 
at t = 10 are practically indistinguishable from an exact 
analytical asymptotic pdf ([7]). The convergence of p(x, t) 
towards /?* (a;) appears to be relatively fast irrespective of 
the chosen /i-driver. 

Although our reasoning is definitely path-wise and all 
data have been extracted from trajectory ensembles, it 
is instructive to visualize generic sample paths. That 
is accomplished in Fig. 3, basically to indicate their 
(paths) qualitative typicalities. The structural impact 
of larger against smaller jumps can be visually compared 
and has been found to conform with standard simulations 
of Levy stable sample paths (with no forces or potentials 
involved), c.f. [27} . 

B. Logarithmic confinement 

1. Quadratic Cauchy target 



-z 2 /2-xz 



-dz 



1 — xz 



-dz 



£12 - £1 
2£i£i2 



[2+ (a; 2 - l)£i£i 2 ] - a; In 



£12 



£1 



(9) 



In the interval z € [£12, £2] we evaluate the integral nu- 
merically. The proposed hybrid procedure (integrating 
analytically in the " most dangerous" small z interval and 
numerically otherwise) permits to speed up the calcula- 
tion drastically and has actually been used in our simu- 
lations. 

The results of our numerical simulations are reported 
in Figs. 1 through 3. We note, that on Fig. 1 the second 
moment oscillates near its equilibrium value 1/2. The os- 
cillations are smoothed out with the growth of the num- 
ber of random trajectories that contribute to the statis- 
tics. A numerical convergence to < X 2 >=l/2 is con- 
sistent with an analytic equilibrium value of the second 
moment of the chosen p*(x) (JT)). The rate of this con- 
vergence is higher for larger /1 <G (0, 2). Clearly, for small 
p, the big jumps arc frequent which enlarges the inferred 
time intervals At in the Gillespie's algorithm, see the tra- 
jectories on left panel of Fig. 3. Thus, the relaxation to 
equilibrium is slow. It gets faster for larger p, when big 
jumps are rare and time intervals At are generically very 
small. 

Fig. 2 displays a probability density evolution, inferred 
from the ensemble statistics of 75000 trajectories. All of 
them have started form the same point x = 0. Although 
the data fidelity grows with the number of contributing 



Let us consider a long-tailed asymptotic pdf which is 
a special a = 2 case of the one-parameter a-family of 
equilibrium (Boltzmann-type) states, associated with a 
logarithmic potential $(x) = aln(l + x 2 ), a > 1/2, see 



p*(x) 



1 



7r(l + a; 2 ) 2 ' 



(10) 



The transition rate ((2]) w$(z + x\x) for any p £ (0,2) 
takes the form 



w<j,{z + x\x) 



C u 



1 



\z\ 1 +»l + {z + xf 



[11) 



Similar to the previous Gaussian case, the simulations 
can be speed up by analytical evaluation of some inte- 
grals. Such acceleration of numerical routines permits 
to handle (in the same timescale) trajectories for much 
longer running times {t ~ 400) than in the previous har- 
monic case (t ~ 10). Each p -driver case (p = 0.5, 1, 1.5) 
will be addressed separately. 

Case of p = 1/2. For z > we need to evaluate 



h/2{x,z)= I _^ dz. 



As 



z 3 / 2 1 + (z + x) 2 



2x + z 



1 1 + x 2 1 

1 + ( Z + X) 2 'z^J 2 " yfl [1 + (X + Z) 2 } ' 



(13) 
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FIG. 1: Gaussian target: Time evolution of the pdf p(x,t) second moment for 25 000 (left panel), 50 000 (middle panel) and 
75 000 (right panel) trajectories. Insets visualize the oscillations smoothing in the asymptotic regime for 10 < t < 15; figures 
near curves correspond to fi values. 




FIG. 2: Gaussian target: Time evolution of p(x,t) inferred from 75 000 trajectories: /i = 0.5 (left panel), n — 1 (middle panel) 
and n = 1.5 (right panel). All trajectories originate from x = 0, i.e. refer to the <5(a;)-type initial probability distribution. 




FIG. 3: Gaussian target. Qualitative typicalities of sample paths for /i = 0.5 (left panel), /i = 1 (middle panel) and n = 1.5 
(right panel). All trajectories originate form x — 0. Right panel comprises two trajectories instead of three (for better visibility). 
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fi/2{x, z) can be written as 

hlM =-^-j M *li; +z?) d , (i4, 

For small z the dominant contribution comes from the 
first term so that the second term can efficiently be eval- 
uated numerically. We note here that although the in- 



tegration in the second term of (fT4|) can be performed 
analytically, the result appears to be quite cumbersome. 
Therefore, for our purposes it is more profitable to inte- 
grate this term numerically. For z < the corresponding 
integral can be expressed as —C^j\i%{—x, —z). 

Case of fi = 1. Here, we need to evaluate the integral 



fi(x,z) 



l + (z + x) 2 z 



dz 

~~2 



In (1 + (z + x) 2 ) + 
I 



■ arctan(z + x) 



2.x- 



InUI 



(15) 



We pay attention to the fact that the integrand is ratio- 
nal, hence the integral can be evaluated analytically in a 
simple form. This permits to use the analytical answer 
(|15[) in our simulations without the need to divide the 
integration range into "small" and "large" z domains. 

Case of fi = 3/2. In this case, we have for z > 



Since 



1 l + x A 



1 / 1 2x + z 



z 5 / 2 l + (z + x) 2 z 3 / 2 \z l + (x + z) 2 
there holds 



(17) 



h/2{x,z) - -^+ (1+a;2)z l/ 2 + 

2xz + 3x 2 - 1 
(I + x 2 )(l + {x + z) 2 ^ 1 / 2 Z ' [ ' 



Again, the third term in ([18]) can efficiently be evaluated 
numerically. For z < we encounter —f3/2(~ x , —z). 

Simulation results arc displayed in Figs. 4 and 5. If 
we compare Fig. 4 with Fig. 1 we see the existence 
of small oscillations in the asymptotic regime about the 
value 1/2. Those from Fig.l were relatively small and 
were quickly smoothed out with the growth of the number 
of trajectories used to extract statistical data. In Fig. 4 
the oscillations are more noticeable and persist even for 
200000 trajectories and more. This is related to much 



slower decay of transition rates (|11[) (determined by slow- 
decaying asymptotic pdf (1101) ) as compared to those for 
Gaussian case ([8]). 

The second moment of the present /9*(x), (|10p . equals 
1 and the convergence towards this value is clearly seen 
in Fig. 3. This convergence is much slower than in the 
Gaussian (harmonic confinement) case which is not a sur- 
prise: (0) and (fTT|)) indicate that the present rate of con- 
vergence should be logarithmically slower. Fig. 5, quite 
alike Fig. 2, convincingly demonstrates a convergence 
of p(x,t) to the asymptotic p*(x). For definitely large 
times around t = 400, p(x,t) and p*(x) become practi- 
cally indistinguishable. Similarly to the Gaussian case, 
the rate of convergence becomes larger with the growth 

of n e (0,2). 



2. Cauchy target 



Now we consider an asymptotic pdf of the form : 



I \ 1 1 



7T 1 + X 

In this case, the transition rate from x to x + z reads 



(19) 



M z + x\x) = J ^ ] j i ^ z ^ x)2 . (20) 

We consider Cauchy driver corresponding to p = 1. The 
transition rate integral can be evaluated analytically. For 
z > we have 



f(x, z) = 




-Vl + x 2 ^J\ + (x + z) 2 + xz In ((1 + x 2 + xz + \/l + x 2 ^J\ + (x + z) 2 ) / z) 



(1 + x 2 )z 



(211 



For z < the outcome is —/(—&, —2). 



In Fig. 6, we report the time evolution of the statis- 
tically inferred p(x,t), its half- width (as second moment 
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FIG. 4: Quadratic Cauchy target: Time evolution of the pdf p(x,t) second moment for 50 000 (upper left panel), 100 000 
(upper right panel), 150 000 (lower left panel) and 200 000 (lower right panel) trajectories. 



does not exist for asymptotic pdf (|19l) ) and simulated 
cumulative probability distributions (CPD) for different 
time instants. An approach to the asymptotic pdf (|19l) is 
clearly seen, together with a convergence of a half-width 
to its asymptotic value 1 . The same convergence pattern 
is observed for CPD which approaches the asymptotic 
function F(x) = \ + 

Comment 3: Displayed empirical (numerically re- 
trieved) curves in Fig. 6 are hampered by certain er- 
rors. The figures have been read from a histogram of 
randomly sampled data. Its partitioning into subinter- 
vals is a source of inaccuracies. In case of a small number 
of intervals, the read-out error would be large, with a size 
of about half-interval length. A finer partitioning (large 
number of small subintcrvals) would still produce an er- 
ror which is close to the half-maximum of the curve. The 
error bound would be smaller or equal to the half-length 
of subintervals corresponding to roughly the same his- 
togram values. One more inaccuracy source in the finer 
partition case comes from the maximum read-out impre- 
cision. Namely, we can have a conspicuous peak, whose 
close vicinity displays much (half or less) smaller values. 
Therefore the partitioning finesse must be slightly opti- 



mized. 



C. Locally periodic confinement 

To set firm grounds for future research it is instructive 
to study our model for more complicated forms of con- 
fining potentials. In view of their physical relevance, it 
is appealing to address an issue of confining (trapping) 
environments with a periodic spatial structure. Here, we 
encounter a major difficulty with a i 1 (i?) integrability of 
the Boltzmann-type weighting function exp(— $). Peri- 
odicity and integrability can here be reconciled either on 
compact sets or by means of locally periodic potentials 
that take a definite confining form (harmonic or poly- 
nomial) for larger values of x € R. Let us consider the 
following asymptotic pdf 

P.(*) = { i _(^_ 4) ' ?, ( 22 ) 
t y x \x\ > 2, 
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FIG. 5: Quadratic Cauchy target: Time evolution of p(x,t) inferred from 200000 trajectories for p = 0.5 (left panel), p — 1 
(middle panel) and p — 1.5 (right panel). All trajectories are started from x — 0. Note scale differences on vertical axes. 




FIG. 6: Cauchy target: Time evolution of pdf p(x,t) (left panel), half-width (HW) of p(x,t) (middle panel). Right panel 
reports the cumulative probability distributions (CPD) for different time instants. Here p — 1 and all data are inferred from 
200000 trajectories, starting from x = 0. 



where C = 3.032818 is a normalization constant. The 
transition rate from x to x + z reads 

M* + x \ x ) = j^rk gx p - ^ x + z ))/ 2 !' ( 23 ) 

where the potential <f> has the form 

We consider fx = 1. To optimize the simulation, here 
we use the same trick of isolating of "most danger- 



ous" small z terms in the integrals involved in the 
Gillespie algorithm. For small z we expand the term 
cxp[(sin 2 (27rx) — sin (27r(a; + z)))/2] in Taylor series. We 
choose e 12 = 0.05. In the vicinity of \x\ = 2 due atten- 
tion must be paid to the proper power series truncation, 
to correctly choose the intervals where integration should 
be performed numerically. For example at x € (1.95,2) 
we have 



The numerators of integrand fractions have been ex- 
panded into Taylor series and (safely) truncated at the 
quadratic terms. 
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FIG. 7: Time evolution of p(x,t) inferred from 200000 trajec- 
tories at fi — 1. The data for 100000 and 300000 trajectories 
(not displayed) do not show qualitative differences. 

Time evolution of the inferred pdf p(x, t) is reported in 
Fig. 7. All sample trajectories were started form x = 
which corresponds to the S(x)-type initial distribution. 
The probability density spreads out with time in confor- 
mity with the trapping (confining) properties of the lo- 
cally periodic enclosure (environment or "potential land- 
scape"). For large running times t=400 the trajectory 
statistics produces data that are indistinguishable from 
those for the asymptotic pdf. We have checked that be- 
ginning from about 100 000 trajectories, further accumu- 
lation of the trajectories number like e.g. 200 000 (dis- 
played) and 300 000 (not displayed) for the data statis- 
tics is inessential. In such cases the curves are almost the 
same, we merely improve a fidelity of the statistics. 



IV. CONCLUSIONS 

If a random process does not admit the description in 
terms of a stochastic differential equation (e.g. Langevin 
modeling), its direct numerical simulation becomes im- 
possible by means of existing popular algorithms. In the 



present paper, for the first time in the literature, we pro- 
pose a working method to generate stochastic trajecto- 
ries (sample paths) of a random jump-type process with- 
out resorting to any explicit (or numerical) solution of 
a stochastic differential equation. To this end we have 
modified the Gillespie algorithm [2l], [22[, normally de- 
vised for sample paths generation if the transition rates 
refer to a finite number of states of a system. 

The essence of our modification is that we take into ac- 
count the continuum of possible transition rates, thereby 
changing the finite sums in the original Gillespie algo- 
rithm into integrals. The corresponding procedures for 
stochastic trajectories generation has been changed ac- 
cordingly. In other words, here we "extract" the back- 
ground sample paths of a jump process, whose pdf obeys 
the transport equation (generalized Fokkcr-Planck dy- 
namics) ([1]), @. We emphasize once more here, that we 
have focused on those background jump-type processes 
that cannot be modeled by any stochastic differential 
equation of the Langevin type. 

Although heavy-tailed Levy stable drivers were in- 
volved in the present considerations, we have clearly con- 
firmed that an enormous variety of stationary target dis- 
tributions is dynamically accessible in each particular 
[i € (0, 2) case. That comprises not only a standard 
Gaussian pdf, casually discussed in relation to the Brow- 
nian motion (e.g. the Wiener process). Among heavy- 
tailed distributions, we have paid attention to the Cauchy 
pdf which can stand for an asymptotic target for any 
\i 7^ 1 driver, provided a steering environment is properly 
devised. In turn, the Cauchy driver in a proper environ- 
ment may lead to an asymptotic pdf with a finite (in fact 
arbitrarily large) number of moments, the Gaussian case 
being included f[l3|). 

An example of the locally periodic environment has 
been considered as a toy model for more realistic physical 
systems. Our major hunch are strongly inhomogeneous 
" potential landscapes" , [l2j , being sufficiently smooth to 
avoid a direct reference to random potentials, 0. Even 
if various mean field data are available in such (exper- 
imentally realizable) systems, it is of interest to have 
some knowledge about the microscopic dynamics (ran- 
dom paths) for the system under consideration. The de- 
tailed analysis of sample path data (ergodicity, mixing or 
lack of those properties) deserve a separate analysis. 

We mention possible generalizations of our method to 
the Brownian motor concept (see, e.g., Rcf. [28| for re- 
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cent review) to include a non-Gaussian jumping compo- 
nent. In those systems it is the properly tailored peri- 
odic "potential landscape" which enforces a conversion 
of a homogenous stochastic process (Brownian motion 
for reference) into the directed motion of particles at 
nanometer scales. That is closely related to the prob- 
lem of so-called sorting in periodic potentials [2!| . Other 
problem to be addressed concerns ultracold atoms in op- 
tical lattices subject to random potentials [3(|, which 
might promising not only from a purely scientific point 
of view, but also with prospects for many technological 
applications. We note that the theoretical description 
of the above mentioned topics relies essentially on the 
Langevin-like equation input. 

Our approach offers an immediate generalization for 
generalizations, where systems with non-Langevin re- 
sponse to external potentials may come into considera- 
tion, along with more traditional ones. What we actually 



need to implement our version of Gillespie's algorithm is 
the knowledge of jump transition rates of those random 
systems only. 

A preliminary work (in progress) shows that an exten- 
sion of our algorithm to higher dimensions is operational. 
In particular, the planar case is worth exploration, pos- 
sibly with more complex "potential landscapes". While 
departing from final comments of Ref . [l2[ we expect that 
the presented methodology can be effectively adopted to 
constru ct op timal random search routines, see in this con- 
nection [3l| . 
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